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We present an almost fully analytical technique for computing Casimir interactions between pe- 
riodic lamellar gratings based on a modal approach. Our method improves on previous work on 
Casimir modal approaches for nanostructures [TJ by using the exact form of the eigenvectors of such 
structures, and computing eigenvalues by solving numerically a simple transcendental equation. In 
some cases eigenvalues can be solved for exactly, such as the zero frequency limit of gratings modeled 
by a Drude permittivity. Our technique also allows us to predict analytically the behavior of the 
Casimir interaction in limiting cases, such as the large separation asymptotics. The method can 
be generalized to more complex grating structures, and may provide a deeper understanding of the 
geometry-composition-temperature interplay in Casimir forces between nanostructures. 

PACS numbers: 



I. INTRODUCTION 

Geometry, material composition, and temperature can 
strongly influence the Casimir interaction [5J between ob- 
jects separated by micron and sub-micron gaps. Recent 
theoretical developments have shown how to compute the 
Casimir force between complex structures using a vari- 
ety of methods [SHE]. Among these, we mention tech- 
niques based on the summation of zero-point energies 
[7HS], which are suitable for high symmetry problems; 
the scattering approach which requires the computation 
of the reflection matrices of the scatterers (TUHEI]; and 
full- wave numerical techniques, that compute the force 
from the Maxwell stress tensor 5, 6, 14 . 

In a previous paper [T] a modal approach was proposed 
to calculate finite-temperature Casimir interactions be- 
tween 2D periodically modulated surfaces. This method 
uses the scattering formula for the Casimir free energy 
and computes the reflection amplitudes of the scatterers 
by decomposing the electromagnetic field into their natu- 
ral modes. The modal approach is based on a plane- wave 
expansion of the fields and a Fourier decomposition of the 
spatial-dependent permittivity of the structures, in the 
same way as done in rigorous coupled wave approaches 
(RCWA) in classical photonics [TS]. The modal method 
is limited to periodic structures, such as photonic crys- 
tals and metamaterials. While other more general nu- 
merical scattering techniques exist, the modal expansion 
provides insight into the different (photonic, plasmonic, 
etc) mode contributions to the Casimir force, thereby al- 
lowing to unveil otherwise hidden balances [16j U-2]- 
Other RCWA techniques, not based on modal methods, 
have been also used by the Casimir community to study 
Casimir forces [3] and nanoscale heat transfer [TS] in grat- 
ing structures. 

In [JJ both the eigenmodes and their eigenfrequencies 



were computed numerically by solving a non-self-adjoint 
eigenvalue problem [191 120) . In this paper we improve 
this previous work by developing an almost fully ana- 
lytical modal approach to compute Casimir interactions 
between ID lamellar grating structures, which is a gen- 
eralization to Casimir physics of well- developed methods 
in grating theory [5TJ [55] . The key feature of our method 
is that the eigenmodes of the grating can be solved for 
analytically without any Fourier expansion of the permit- 
tivity, while the eigenfrequencies are solutions to a simple 
transcendental equation. Analytical expressions for the 
eigenfrequencies can be found in some limiting cases, such 
as for perfectly reflecting gratings, and for the low fre- 
quency limit of real material gratings, described by sim- 
ple Drude or plasma permittivities. The quasi-analytical 
modal approach also allows us to exactly demonstrate 
some properties of the scattering operators and to de- 
rive expressions for the Casimir interaction in some limit- 
ing cases, such as the large distance/low frequency limit, 
and the behavior of the force at high temperatures. The 
method can be generalized to more complex structures 
beyond ID lamellar grating, and can also provide a de- 
tailed framework for the analysis of other fluctuation- 
induced interactions in nanostructures, including thermal 
emission and near-field heat transfer. 

The general set up is similar to that of [TJ, which we 
briefly outline here. Within the framework of the scatter- 
ing approach, the calculation of the Casimir free energy 
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is essentially reduced to the calculation of the scatter- 
ing matrices of isolated objects. The symbol Tr indicates 
the trace over spacial and polarization degrees of freedom 
[TJ. Here /3 = is the inverse temperature, X rep- 
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FIG. 1: A schematic representation of the scattering of the 
electromagnetic field on a ID grating. 

resents translation matrices that depend of the distance 
a between the gratings, and 1Z are their reflection matri- 
ces. All these matrices are evaluated at the Matsubara 
imaginary frequencies w; = = i2nlkBT/h |23j . and 
the primed sum indicates that the I = term has half 
weight. The arrows under the reflection and translation 
matrices indicate the direction of propagation of light - 
for example, is the reflection on the left grating for 
light propagating from left to right. The translation ma- 
trices are diagonal in a plane-wave, Rayleigh basis (see 
PQ for explicit expressions). 

Our goal in the rest of the paper is to compute the 
reflection matrix of an isolated ID lamellar grating with 



the quasi-analytical modal technique. In the following 
we will analyze the scattering properties of a ID lamellar 
grating of depth d and period p = p\ + P2 , where pi is 
the width of the grooves and P2 the with of the teeth 
(see FiglTl). We divide the ID lamellar grating into three 
regions: (i) the homogeneous, vacuum region above the 
grating, z > 0; (ii) the region z < —d below the grating, 
filled with homogeneous medium of permittivity e(w) and 
permeability /i(w); and (iii) the grating region — d < z < 
0, where the space is filled with the modulated medium 
e(x;ui) and /i(x;lu) describing the ID lamellar grating. 
In each i-th region (i = v, vacuum region; i = g, grating 
region; and i = m, bulk medium region), the solution of 
Maxwell equation can be written as 



/E x (x,z)\ w 

\Hy(x,z)/ 

= ^A^Y^^A^'V^'*' 2 , (2) 

where we have considered the four independent trans- 
verse field components (the remaining two components 
can be directly calculated from the previous four) . More- 
over, since the system is invariant with respect to trans- 
lations along the y-direction, eigenmodes have a e zkyV 
plane wave form, which we omitted in the above equa- 
tion. Y^[aj,Ai a,<) ] is a column vector describing the x 
dependence of the eigenvector with corresponding eigen- 
value Xu'^, where v labels the eigenvalue and s denotes 
one of the two possible polarizations. A^?'*' are complex 
amplitudes, to be determined by imposing the following 
boundary conditions at the interfaces (continuity of the 
tangental components of E and H at the interfaces) : 



J2 A t g)Yis ' g) [x,K] 
Ai u' m) Y(s,m) [m, K]e- iX " ,m) 



where we have dropped the superscript for the eigenval- 
ues argument of the eigenvectors because they are the 
same as the eigenvectors. Using properties of the eigen- 
vectors (see below) it is possible to derive the Alf in 

terms of the A^ ,m ^ or vice-versa, and then extract the 
scattering operators of the grating. 

The expressions for the eigenvectors and the corre- 
sponding eigenvalues play a central role in our derivation. 
In the following we will focus on the evaluation of these 
quantities for the grating region (i = g). The results will 



= J2Ai^Y^[x,X v ], (3a) 
d = H A l 8 ' 9) yK9) [*> A„]e- iA ^ 9)d , (3b) 
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be also valid for the bulk and vacuum homogeneous re- 
gions. For this one simply needs to take the limit of no 
modulation (e(x;uj) = e(u>) and (j,(x;lj) = /z(w)). 



II. NON-SELF-ADJOINT EIGENVALUE 
PROBLEM 

In this section we give the mathematical details on how 
to solve Maxwell equations in a ID modulated magneto- 
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dielectric region. The modulation is lamellar along the 
x— direction (the electric permittivity and the magnetic 
permeability modulated along that direction), and the 
grating is invariant along the y— direction. For the pur- 
poses of finding the eigenmodes in the grating region, we 
assume that the system is also invariant along the z— 
direction [5TJ [55]. Using the invariance in y, it is possi- 
ble to write Maxwell equations V x E — iw/iH = and 
V x H + iweE = in the following form 



where we already eliminated the ^-component of the elec- 
tric and the magnetic fields. For the sake of simplicity 
we will also be measuring all frequencies as wave vectors, 
so that lu/c — > lo. 

The previous system of equations can be solved by sep- 
aration of variables, by writing 



d z E x 



d z H v 



<>, ( / ' " II . , 

. we 



d x - — d x — iuifi 
iuje 



k 2 k 
d z E y = — - — H x H d x H y , 

dz H x = d x ( — E, 



d x - — d x — iojfi 

lUlfJ, 



k 2 k 

— E- E 

lUJfJ, UJjX 



H y , (4a) 
(4b) 
v, (4c) 
(4d) 



E, 



(E x {x,z) 

Ey(x,Z) 

H x (x,z) 

\H V (X,Z)J 



(E x [x,XY 

Ey[X,X) 

H x [x,X] 

\Hy[x,X}/ 



where Y[x, X] satisfies 



Yfx,Ale 



(5) 
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XY[x,X] = 



i ujfi 
_ *i 



ujfi i 



i UJE 



[d x ±d x +ujfj]\ 



uje i 



Y[x,X], 



(6) 



which is a eigenvalue equation for the eigenvector Y with 
eigenvalues A. Here k 2 = k 2 — k 2 and k 2 = fiecu 2 , and for 
simplicity we have omitted the spatial and frequency de- 
pendency of the permittivity and permeability functions. 

The 4x4 matrix equation ^ can be easily transformed 
into a 2 x 2 second order differential equation either in the 
E- or _ff-components only. As a further simplification we 
decompose the fields in two independent polarizations. 
We will define "e" or "h" polarizations, for which the 
x-component of the electric or magnetic field vanishes 
respectively, i.e. E x — and H x = 0. Using this de- 
composition it is possible to show that the 2x2 matrix 
equation decouples into two one dimensional second or- 
der (in general non-self-adjoint) differential equations for 
the y-components of the fields, namely 



<7 {S \x)d x 



<t( s )(x) 



d x + k 2 {x) 



U (s) [x,X] = X 2 U {s) [x,X], 
(7) 



,h, a^(x) = fj,(x), ^(i) = e(x), 



where s 

E y e \ and = H y h \ The corresponding eigenvalue A 
will therefore also depend on the polarization s. 
Given Eq.([7|), from Maxwell equations we get [55] 

J^)[x,X} = -^—^d x U^[x,X], (8) 



where S^ e 



-1, 5W = 1, J( £ ) = H y e \ and jW = 

E x n '. Solving (B and using the solutions in |s| one can 
therefore find the y components of the eigenvector Y, 
and from them, using again Maxwell equations, the x 
components, given by 
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iX[iuj\ 
iXeui k y d x 



kyO x 



Finally one can write the eigenvectors for the two polar- 
izations 



Y&[x,\] 



Y(*)[i,A] 



( ° \ 

I k y d x U^[x,\] 



ky d X U (H) lx,X] 







(10) 



V U { V[x,X] J 



The matrix differential operator in the r.h.s of Eq.([6j) 
is clearly not Hermitian. Therefore the eigenvalue prob- 
lem (|6| is non-self-adjoint and the eigenvalues A are in 
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general complex [Tnj[2U]. Note that, since the matrix is 
non-symmetric, this remains true even if we consider a 
non dissipative material. The existence of such complex 
values is associated with the presence of evanescent fields 
in the structure [21] • 

Following the theory of non-self-adjoint differential 



equations [T5J [20], one also needs to find the adjoint 
eigenvectors, which are bi-orthogonal to the eigenvectors 
in Eq.( 10 ), in order to completely characterize the math- 
ematical description. Indicating them with Y[x, A], one 
can show that they are solutions to 



AY[af, A] = 



uje* i 



uifj,* i 

d x ^d x +uje* 






LJfl* 

d x by 
i ujfi* 



Y[x,X] 



(11) 



The two eigenvalue problems ([6]) and ( 11 ) have the same 
eigenvalue spectrum [T^l HO]- One can find the adjoint 
eigenvectors employing the same strategy used above. 
One gets 



c V- e) [x 



Y ie) [x,X] 



i/i*(x) 



Al\ 



/Was) 



X 2 +k 



Vv( fi )[x,A] 



Y W [x,A] 



/ 

i-j%VW[x,\]\ 

"o 

k y a J: v w [x,x\ 

A 2 +feJ ie*(x) 
V (h) [a;,A] 
V 6*(*) 



(12) 
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where the function yW[a;, A] satisfies the differential 
equation 



[a^(x)]*d x 
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V (s) [x,A] = A 2 V (s) [x,A], 
(13) 

which is the adjoint equation of The eigenvectors 
and their adjoint are bi-orthogonal, i.e. 



(Y" (s) [A 7 ]|Y( s ')[A y ]> = / dx 

J-p/2 
— ^7,7' $s,s f • 



p/2 



Y (s) [a;,A 7 ]t • Y^[x,\y) 



a( s \x) 



(14) 



In deriving the previous expression we used that the func- 
tions IA and V satisfy themselves the bi-orthogonality re- 
lation 



^^V^[,,Ay')[,,A yj 



1 



p/2 



<t( s )(x) 



7,7' 



(15) 



The choice of the normalization factor of V^ s ^ and K^ s ' 
allows us to have V^*[a;, A] = U^[— x, A], which will be 
particularly handy for the forthcoming evaluations |22j . 



It is interesting to note that, from the symmetry in the 
equations in Q and (11 1, one can also show that the bi- 
orthogonality relation ( 14 ) also has a physical meaning. 



It is directly connected with Poynting vector reciprocity 
theorem, and therefore with the energy flux along the z 
direction. 



III. EIGENMODES AND EIGENVALUES 

The approach described in the previous section re- 
quires the solution of the second order differential equa- 
tion given in Eq. Q. Because of the periodicity 
of our system the functions and their derivatives 
must satisfy pseudo-periodic (Bloch-Floquet) boundary 
conditions, i.e. U^\p/2,X] = e la ° p U^ [-p/2, A] and 
d x U^\p/2,X] = e iaoP d x U^[-p/2,X]. The previous re- 
quirements define a (in general non-self- adjoint) scalar 
eigenvalue problem. In this section we give the form of 
the eigenfunctions [x, A] and the corresponding eigen- 
values A, both for the grating and homogeneous regions. 



A. Eigenfunctions in grating region 

Within the modulated region (—d < z < 0) of the 
ID lamellar grating, the permittivity and permeability 
functions are 




(16a) 
(16b) 



Using these expressions in the differential equation ^ 
one can find its solutions [x, A] must have the follow- 
ing expression 
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WW[X,A] = {cos(a p/2)^ ) [x,A] +i S in(a oP /2)^[ S ,A]} 



(17a) 



where ao describes the x component of the wave vector limited within the first Brilloin zone {—n/p < k x = <xq < tt/p), 

ipi s) [x,X] = c/> { e s) [x,\}/<t> ( e ) [p/2,\} and ^[as.A] = [a:, A]/<^ s) [p/2, A]. The functions (f> { e } [x,X} and ^[x.A], even 
and odd in the variable x respectively, are given by 



,A1 = 



cos( 7 ia:), for \x\ < ^ 
cos ( 7l f) cos ( 72 [|x|-f]) 

s mfax) 



^72 



sin( 7l f)sin( 72 [H-f]), for f < |z| < § 



4 S) M] 



7l 



for |z| < f 



si„( 7l ^) co S ( 72 [b|- ^]) + ^r^ co S ( 7l sin[ 72 (|x|-^)] 

sign(x) ^ 



(17b) 
(17c) 



for f < \x\ < | 



where 7 f = ei/iiw 2 — (fc 2 + A 2 ) and 7 | = C2/i2W 2 — (A; 2 + A 2 ). The normalization constant C^ S '(X) is given by 



C (s) (A) = 



cos 2 (a p/ 2 ) 



da; 



^\x,Xf 
a( s )(x) 



sin 2 (a p/2) 



d.r 



yl s) [x,A] s 
a^(x) 



(17d) 



We note that all the above functions, being even in 71,2, 
do not depend on the definition (sign) of the square 
root, i.e., do not contain any branch cut. We also note 
that U^[x,X] = U^[x,—X]. The previous expressions 
( 17) are fully determined only once the eigenvalues A are 



known. 

B. Eigenvalues in the grating region- General 
properties 

Imposing the pseudo-periodic boundary conditions (on 
the function and its derivative) it is possible to show that 
the eigenvalues are the solution of the following transcen- 
dental equation [52] 

= L>( s) (A) = - cos(a p) + cos(pi 7 i) cos(p 27 2) 

-T«(A) sin(p l7l ) sin(p 272 ), (18) 

where 



T (S \X) 




(19) 



Eq ( 18 1 clearly shows that the eigenvalues depend on 



the frequency u, and the two wave- vectors ao and k y . 
We can deduce the following properties, valid for both 
polarizations (we will drop the superscript s): 

• D(X) is quadratic in A and, therefore, if A is a so- 
lution then —A is also a solution. 

• D(X) is even in 7, which implies that the solution 
is not affected by the sign of the square root. 



• -D(A) is even in ao and k y , which implies that A is 
an even function of the two variables. 

• For complex frequencies lj = (, A(C) = — A*(— £*), 
which implies that A is a pure imaginary quantity 
for imaginary frequencies u> = i£. 

• At high frequencies the ultraviolet transparency of 
all materials (e = fj, = 1 for uj — > 00) implies that 



X(lo — > 00) = ± 



\ 



h 2 

Ky 



a 



2-kv 
P 



(20) 



with veZ (see also Section V A | 



• The solutions of -D(A) = form an infinite, numer- 
able set of complex numbers. 



C. Special case: homogeneous media 

The previous formalism also applies to the homoge- 
neous region of space. This particular limit is reached by 
imposing o\ = er 2 = a in the previous expressions, which 
implies 7 i = 7 2 = 7 . The corresponding transcendental 
equation becomes then 



D(X) — cos( 7 p) — cos(aop) = 0. 



(21) 



There are two possible sets of solutions for 7, namely 
7 = ±a + = ac Vi ±, where v 6 Z. The corresponding 
eigenvalues are 



X v = ±^Jefiuj 2 - (k 2 y + a 2 j± ), 



(22) 



G 



and are the same for both polarizations s — e and s = h. 
From a comparison with the usual grating theory each 
value of v corresponds to a specific Rayleigh order. How- 
ever, it is important to note that for the eigenvalues A 
the set of solutions with the + sign gives an identical re- 
sult as the set of solutions with the — sign. This means 
that, for v = 0,±1,±2,..., one needs to consider either 
one or the other. Instead, both set of solutions must be 
considered if we limit v = 0, 1, 2, . . . (the v = solution 
must be counted only once). The eigenfunctions are 



UW[x,\ v ] = (-l)«y ^-e^l+^h (23) 
which are the usual plane wave Rayleigh modes. 

IV. SCATTERING OPERATORS 

In this section we describe how to compute the re- 
flection matrices of the nanostructure within the modal 



approach. What follows is similar to what is presented 
in [1] with the important difference that now we have 
analytical expressions for the eigenvectors. We introduce 
a transfer matrix that relates the field amplitudes of the 
vacuum and bulk regions, and express the scattering op- 
erators of the grating in terms of the transfer matrix. 



A. Transfer matrix of the grating 

We start by splitting the complex eigenvalues A into 
two subsets: (a) eigenvalues with positive real part, with 
corresponding eigenvectors called "right" eigenvectors, 
and (b) eigenvalues with negative real part, with corre- 
sponding "left" eigenvectors. Each subset is then ordered 
by the increasing moduli of the eigenvalues, the smallest 
eigenvalue denoted as A„ = o for the subset (a) [— X u =o for 
subset (b)], X u= i [— A„=i] for the next eigenvalues, etc. 
Although this ordering is not unique, the end results are 
not affected by our choice. 

We define the fundamental 4 x oo matrices 



3,W = (| Y W[a£2]>, |YW[A^]), |Y<«0[a£>]>, |YW[A£l?]>, 



W 



|YW[-Al e J>, |yW[-a££>]>, |Y(e)[-A(fi]>, |YC0[-A£$]>, • • •) , 

I 



(24) 



where we dropped the x dependency of Y^[a;,A]. Wc 
recall that the index i indicates the region under consid- 
eration (i — v,g,m). We also define the oo x 4 matrices 
formed by the adjoint eigenvectors 



as 



F«(x,z) = (£ w ,3> w ) 



'•pw 

V (l } 



(27) 



where is a column vector formed by the amplitudes 
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Wt 



((Y {e) m}\\ 

(y (h) [x^}\ 

(Y (e) [A^J]| 
<Y ( 'V»>]| 

J 



Wt 



f(Y ie) 




(Y {h) 


-a£^]| 


(Y (e) 


-A^]| 




-A^?]l 



V 



A<£ A in" 



By applying the boundary conditions ([3]), one gets the 
following relation between the amplitude coefficients 



/ 

(25) 
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(28) 



Using the scalar product defined in (pk we have y Wt . where we have defined the vectors as the field am 



3£W = I and £«t ■ V_ W = 0. Let us define also the 



*7 ^. ^ 

diagonal propagation matrices 



V (l) (z) = diag[e iA -° 



z J\ y 



plitudes in the bulk media multiplied by the correspond- 
ing phase factors e~ lX » ' >d , i.e. the field amplitudes at 
the bulk/grating interface (z = —d). The grating trans- 
fer matrix O is a 2 x 2 block matrix, with each of the 
blocks defined as 



£ {l) (z) =diag[e 



-iA 



(e,0 
= 2 



-iA 



(26) 



which clearly verify [P^\z)} 1 = J^_{z) = T^{—z). Using 



the previous definitions, the field in eq.([2]) can be written 



On 


_ -y (m)f 


G(d) 


1 {V \ 


(29a) 


Ol 2 


_ -y (m)f 


G(d) 




(29b) 


Q21 


— V Mt 
4 


G(d) 




(29c) 


e 22 


_ -y(m)f 


&{d) 




(29d) 
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where the "grating" operator G(d) is 



G(d) = ^ (9) • £_ {9) (-d) ■ £<»>t + y(g) . ?>0)(-d) . y Grtt 
= £|Y«[AH])<Y (s ^)]| 



-iA, ( ,"' B) d 



+ (A^- 9) -> -Xi s - 9) ) = G (e) (d) + 



(30) 



Thus, the operator G(d) is a decomposition in function 
of the polarization as well as the right and left eigenvec- 
tors. The grating operator describes the propagation of 
the electromagnetic field through the grating, and it is di- 
rectly related to the Green tensor of the electromagnetic 
field in the modulated region. 

Once we have obtained the theta-matrices we can get 
immediately the scattering operators: 

g = -e 22 1 -e 21 , (3ia) 
^ = e 12 -e 22 1 , (3ib) 
j_ = 6u - e u • e^ 1 • e 21) (3ic) 

% = ©22 1 - (31d) 

The logic behind the previous expressions is simple: 7£ 
is the matrix that gives the field amplitude A^- v > in terms 
of AW. Similarly T connects £S m ^ with AW, etc. As 
one can see, the derivation reflection and transmission 
operators requires the inversion of the matrix 022 , which 
is called the pivotal matrix 22J. This matrix contains 
important information about the scattering properties of 



the grating. Indeed, the zeros of its determinant are con- 
nected with the resonances of the scattering operators. 
A simple check of this property can be found in the next 
section, where the resonances are the surface plasmons 
for a plane metal-dielectric interface. Despite being for- 
mally simple, the inversion of this matrix may pose nu- 
merical problems because the matrix is sparse, can be 
singular (ill conditioned), and may lead to numerical in- 
stabilities. We will see in the last section how one can 
skirt this problem when computing the Casimir energy. 
As a last remark let us notice that, from the properties of 
the eigenvalues and of the eigenfunctions, it follows im- 
mediately that all scattering operators are symmetric in 
«o and k y . This information will simplify the calculation 
of the Casimir interaction at the end of this paper. 



B. Special case: planar interface 

To validate the previous approach and clarify how the 
actual calculation works, let us consider the simple case of 
a planar interface between two homogeneous media ( "to" 
and 'V'). In this case we should recover the expression 
for the Fresnel reflection amplitudes in the e and h po- 
larization basis. For d — the operator G becomes the 
identity, simplifying the expression of the theta-matrices 
(291, which become block matrices, with each block hav- 
ing a dimension 2x2. Using the expressions for the 



eigenvalues (22 1 and eigenfunctions (23) in the homoge- 



neous regions, we obtain all the elements of the pivotal 
matrix: 



(ee) _ ( i i \^ [\P? + ki . 



[9 



[6 



2 V^™ 3 ^ [Ai m) ] 2 + fc,V 
Ve(^W / 1 1 \j m) [\i v) ] 2 + k 2 z \ 

{eh) _ yV^eW k z a v 



22 \iv 



2 



k z ot. v 
e M At («) w A« 




As before, a v = cto + 2ixv /p. The other matrices can be immediately derived from the previous expression by 
accordingly changing the sign of A. For example for QJ 2 , X m — > — X m while for Q 2 i > X v — > etc.. This also means 
that all O-matrices are block diagonal with each block being 2x2 and, therefore, the same occurs for the reflection 
operator. In the special case k y = the (e, h) polarization basis coincides with the usual transverse electric (TE) 
and transverse magnetic (TM) polarization basis. In this case the blocks and therefore the reflection operators are 
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diagonal and, as an example, we have 



V 



fj, v A^+M™ A£ 

o 



{"A"-e'°A; 
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where r TE ' TM are the usual Fresnel reflection amplitudes. 



(33) 



V. EIGENVALUES: ANALYTICS AND 
NUMERICS 

It should be clear from the previous section that the 
key element to calculate the scattering operators are 
the solutions of the transcendental equation. In or- 
der to study them both analytically and numerically, 
it is convenient to define the variable r\ = 7 2 , write 
72 = V + [M2( w ) e 2(w) — ^i(a;)ei(a;)]w 2 , and re-write the 
transcendental equation in terms of the variable 77 as 
D(rj) — 0. The advantage of doing this is that, in contrast 
to ( 18 ), this new equation does not depend on k y , thereby 



reducing the dimensionality of the space where the so- 
lutions are defined. Once we solve for 77, we obtain the 
original eigenvalues A using A 2 = ^ 1 (w)e 1 (w)w 2 — (ky+rj). 

In general, the solutions of the transcendental equa- 
tion D(rj) = must be searched for numerically. This 
task is complicated by the fact that, for real physical fre- 
quencies, they are complex numbers. However, since the 
Casimir free energy ([I]) is given as a sum over the pure 
imaginary Matsubara frequencies, in this section we con- 
sider the solutions of the transcendental equation already 
at imaginary frequencies, namely 



= L> (s) (?7) = - cos(a p) + cos^iy^) cos(p 2 \A? 

1 ( Vv - R<0 - 1]£ 2 , 



m) - 



im(pt y/rj) sm(p 2 \J-q- [e(i£) - 1]£ 2 ) , (34) 



where, for simplicity, hereafter we specialize to the case 
where one of the medium is vacuum (ei, (i± = 1) and the 
other has no magnetic activity (e^ = e, f/,2 = 1)- Our 
derivations and discussions below can be generalized to 
other grating configurations, where, for example, instead 
of vacuum we consider other materials, such as dielectrics 
or semiconductors. We also recall that in the previous 



expression, <J^\i£) = 1 and cr^'i^) = e (*£)- One can 
analytically show that on the imaginary frequency axis 
u> = i£, the solutions for 77 = rj[£, eta] are non-negative, 
real numbers. The eigenvalues are the n p urely imagi- 
nary quanties, A = ±i» /£ 2 + fc 2 + rj. Eq.(34l is the main 



equation in this work, that we shall study in 



detail below. 



Drude (en) and plasma (e p ) permittivities: 



(35) 



where uj p is the plasma frequency and 7 the dissipation 
rate. In the following we study the high and low fre- 
quency behavior of the eigenvalues for both permittivity 
models. 

As we stated in Section II, for £ 3> oj p the ultravio- 
let transparency of metals implies that the Drude and 
plasma models share the same set of eigenvalues, inde- 
pendent of polarization. The large eigenvalues (rj 3> w 2 ) 
have the form 



A. Drude and plasma models for metallic gratings 

Depending on the range of frequency, the dielectric 
model, and the polarization it is possible to find ap- 
proximate analytical expressions for the eigenvalues in 
some limiting cases. For simplicity, we will consider here 
only two model dielectric functions of metals, namely the 



? ? (£ -> 00) = (a + 2irv/p) 2 , 



(36) 



with i/€Z, while for rj < w 2 their values must be found 
numerically. 

At low frequencies the eigenvalues depend on polar- 
ization, and they are different for the Drude and plasma 
models. We will call "low frequency" different regions for 
each of these models: for the Drude mode it corresponds 
to £ <C 7, while for the plasma model to £ <C u> p - In the 
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region 7 <C £ ^ w p the solutions for the two polarizations 
behave differently, but the plasma and the Drude model 
give similar expressions. In the region £ <C 7, absent in 
the plasma model, the Drude model describes a regime 
where the electromagnetic field undergoes a diffusive dy- 
namics (25] [26] • We now consider the two polarizations 
separately. 



1. s = h polarization 

In the limit £ <C uj p , assuming that ?7[£,ao] is constant 
or goes to zero slower than [e(z£) — 1]£ 2 , the s = h value 
of the term in the big parentheses in the second line of 
Eq.(34| is much larger than one. Then one has to look 

for solutions of sin^i^/rj) sin(p2\A7 — [ e (*£) — 1]£ 2 ) = 0- 
Two sets of solutions are possible: the first 



pi 



(37) 



(y G Z, and v 7^ 0) does not depend on the permittivity 
model and describes modes vibrating within the grooves 
(see fig. [3]); the second one is given by 




«+7 



(plasma) 
(Drude), 



(38) 



and describes modes vibrating inside the teeth (see fig. 
[3]). The difference between the two dielectric models is 
evident in the limit ( « 7. For the plasma model all 
the solutions are always distinct. On the contrary, for 
the Drude model degeneracies are possible: for certain 
frequencies £ there are values of v that make the eigen- 
values of Eq. (l37l) identical to the ones of Eq. ( 38 ) , and in 



this case an alternative approach must be used to search 
for the solutions (see the end of this Section) . 

For 77[£,a ] going to zero faster than [e(i£) — 1]£ 2 for 
£ — >• 0, one can no longer neglect the terms in the first 
line Eq.(34). In this case, for £ <C lo p in the plasma 
model one can approximate r\ — [e(i£) — 1]£ 2 ~ — uj p , and 
expand up to the second order in 77 the terms cos^iy^) 
and sin(piy^). Solving the resulting equation one gets 



for the smallest eigenvalue 



(h) 

7 ?i/=0,pla 



2£ 



, cosh(p2^p) — cos(aop) 
u p pi sinh(p 2 w p ) 



(39) 



which describes a mode resulting from the coupling of 
surface plasmons living on the walls of the grooves. For 
the Drude model 77— [e(i£) — 1]£ 2 goes also to zero for £ — > 
0. Expanding to the second order in rj the corresponding 
trigonometric functions, and solving for rj one gets 



(40) 



(h) f [1 ~ cos(a p)] 1 p 2 £ 

Hence, t/q p lasma goes quadratically to zero with the 

frequency, while the corresponding power law % Drude 
strongly depends on the value of a®. This last feature 
will be relevant in the numerical evaluations below, in 
particular in the calculation of the zero frequency limit 
of the reflection operators. 



2. s = e polarization 

Let us consider now the low frequency behavior of the 
eigenvalues in the case of e-polarization. For the plasma 
model one can see that Eq.(34) does no longer depend 



on the frequency, and in consequence the correspond- 
ing eigenvalues are frequency-independent and coincide 
with their high-frequency limit. The eigenvalues must be 
found numerically, the large ones being approximately 
equal to (36). For the Drude model Eq.(34) becomes 



identical to the one for vacuum. The solutions are then 



Vd 



Drude « 7) - 



(±a Q 



2ttv 
P 



fovv = 
for v ^ 0. 



(41) 



In this case degeneracies happen at the center (olq = 0) 
and at the border (ccq = k/p) of the Brillouin zone [S7J 



(see Section VII for the impact of degeneracies on the 
calculation of the Casimir interaction). Expanding the 
transcendental equation ( 34 ) to second order in rj around 



the solutions (41 1, and solving for rj one gets 



(<0 _ 

^Drudc — < 



d r ,bM( v )-j\d n b^)(T } )]*-ib(e)( v )d*bi*){' n ) 
(±«o + ^ 2 



for v = 



for v^0. 



(42) 



J?=(±a +2iry/p) 2 



It is possible to show that in the limit ao — >• 

(e) _ 2 , P2 £ 

7 ^=0,Drudc ~ a + — -■ 



(43) 



B. Numerical solution for eigenvalues 

We solved numerically the transcendental equation 
( 34 ) using Mathematica, employing as seeds for the roots 
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FIG. 2: Numerical solution of the transcendental equation 
( 34 1 for the Drude model. Only the lowest eleven eigenvalues 
are shown, (a) s = h polarization. The dotted line corre- 
sponds to the smallest eigenvalue rj^2 . Solid lines are the 
eigenvalues obtained using as seeds the expression ( 37 1 [u go- 
ing from 1 (bottom curve) to 4 (top curve)]. Dashed li nes are 
the eigenvalues obtained using as seeds the expression) 38 1 [u 
going from 1 (bottom curve) to 2 (top curve)], (b) s — e po- 
larization. The dotted line corresponds to the smallest eigen- 
value T]„2 . Solid and dashed lines are the eigenvalues ob- 



tained using as seeds the expressions ( 41 \ for the two possible 
signs [v going from 1 (bottom curve) to 4 (top curve)]. Dashed 
lines are the eigenvalues obtained using as seeds the expres- 
sion (38 1 [v going from 1 (bottom curve) to 3 (top curve)]. 
Parameters are pi — 160 nm, P2 = 90 nm, and ao = 0.5n/p. 
The optical parameters chosen for these plots are uj p — 8.39 
eV, 7 = 0.043 eV; the general structure of the curves remains 
unchanged for other choices of Drude parameters. 



of D^ 8 '(rj) — the asymptotic analytical expressions for 
the eigenvalues described above. Hereafter we will fo- 
cus on the results for the Drude model, postponing the 
results for the plasma model and their comparison for 
future work. In figure [2] we show the numerical solutions 
for specific values of the geometrical parameters of the 
grating. As seen in the figure, all s — e eigenvalues bend 
down at low frequencies, while only the s — h eigenvalues 
obtained from the seeds ( 38 ) show the same trend. This 
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FIG. 3: Spatial structure of the electromagnetic modes in the 
grating region in units of the plasma wavelength A p = 2tt/u} p , 
for k y — 0, ao = 0.2-7r/p for the second Matsubara frequency 
£ = AnksT/h atT = 300 K. The curves represent the modes 
with v = 2, corresponding to the two seeds in Eqs. (37 1 and 
( |38[ ), and to the two seeds (± solutions) in ( |41[ ). Our choice 
k y = implies that the e and h polarizations decouple, and 
that E x and H y depend only on the /i-polarization, while E y 
and H x only on the e-polarization. For the /i-polarization two 
categories of modes exist: the first mainly vibrate within the 
grooves, and the second mainly within the teeth (this is par- 
ticularly clear for H y ). In agreement with Maxwell equations, 
the component of the electric field along the modulation direc- 
tion is discontinuous at the groove walls while the remaining 
ones are all continuous. The numerical values inside the plots 
indicate the effective refractive index at imaginary frequency, 
n e g{i£) = fcz/C for the corresponding mode (the value on the 
left corresponds to the full line mode, the one on the right to 
the dashed line mode). The Drude and grating parameters 
are the same as in the previous figure. 



For the h polarization the eigenvalues obtained from the 
seeds ( 37 1 change smoothly and they cross the other set 



of h eigenvalues for some values of £, showing degenera- 
cies. Both sets of eigenvalues of the h— polarization are 
almost insensitive to the value of ao, while this param- 
eter becomes relevant at large imaginary frequency. In 
contrast, the eigenvalues of the e-polarization are very 
sensitive to the value of cvo- In n gs- [3] and [4] we plot the 
spatial profile of the eigenmodes corresponding to the 
previously discussed eigenvalues. 



VI. DETAILS ON THE CALCULATION OF THE 
MATRIX ELEMENTS 

Now that we have described the calculation of the 
eigenvalues and the eigenvectors, let us proceed to the 
computation of the theta-matrices (29 1. All these matri- 



ces involved in the calculation of the Casimir free energy 
have a similar form, namely a collection of 2 x 2 blocks 
with elements coupling the two polarizations: 



behavior is due to the dissipative nature of the metal. 



■i)%]\Gw\YV>*)[(-iy\ v ]). 



(44) 
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FIG. 4: Density plot of the intensity of the electromagnetic 
field within the modulated region. Upper box corresponds to 
\E X \ 2 and lower box to |_B H | 2 . In each of the plots the horizon- 
tal axis is the modulation a;— direction and the vertical axis is 
the invariant z— direction. The parameters are the same of 
the previous figure and, therefore, \E\\ depends only on the 
/i-polarization (eqs.( |37| l and ( |38[ ) ) while \Ey\ depends only on 
the e-polarization (eq.1 41 1). The first five modes for each 



polarization and for each set of eigenvalues are represented. 
There is only one zero mode per polarization (see discussion 
in the text). For the parameters chosen here the zero mode 
is almost constant for the /i-polarization. 



From the expression for the grating operator ( 30 1 and of 



the eigenvectors, it follows that one of the key elements 
of our approach is the calculation of the overlap between 
the eigenvectors describing the field in the grating region 
with the eigenvectors characterizing the field in the two 

homogeneous regions, namely (Y S [AS*' m '']|Y^ '[A^ ,9 ' ) ]) 

and (Y W [\l a ' 9) ]\Y ( - a ">[\%' v) ]), and eventually the ex- 
plicit calculation of the following integrals 



UM[x,\] iotx 



dx d x U^[ X ,X] e -^ 



<t( s )(x) 



(45) 



It is interesting to note that since U^[x,X], defined in 
(17a I, is a combination of trigonometric functions oscil- 



lating with frequencies 71 and 72, the above integrals are 
large when a = ±7, (i — 1, 2). Physically speaking this 
relation describes the x-component momentum matching 
between the electromagnetic wave coming from the ho- 
mogeneous regions and the wave propagating in the grat- 
ing region. The above integrals can be done analytically, 
but the resulting expressions are long and cumbersome, 
so we do not report them here. 

It is interesting to consider some special cases. For 
instance, for k y = the scalar products between vectors 
with different polarizations vanish. As a consequence the 
polarizations decouple and one can show that the blocks 
of the theta-matrices become diagonal. Under a transfor- 
mation that generates an even number of permutations 
of rows and columns, the transfer matrix can be written 
in a block diagonal form as 



0(ky = 0) 



0(e 



; ) 

Q(hh) 



(46) 



From Eqs.(31) it immediately follows that all scattering 
operators are 2x2 block diagonal. In the case where the 
Drude model is used to describe the optical properties of 
the metallic grating, some interesting information can be 
obtained for the reflection operator in the limit £ — > 
for the e-polarization. In this limit = U^om because 
the eigenvalues (41) are identical to the ones of vacuum, 
i.e., the grating modes match the ones of the vacuum 
region, and therefore the electromagnetic field effectively 
does not see the grating modulation. The properties of 
the e-polarization allow to directly connect this result to 
the Bohr- van Leeuwen theorem 128] [29]. 



Decomposing the operator G over the two polariza- 
tions, in the limit £ — > and arbitrary k y we can write 



e 7 „(£ - 0) 



(Y (e) [A 



(e,m) 

7 ' 



]\Q(h) | Y W [A^ l,) ]) (Y (e) [A^ e ' m) ]|GW |YW [Al M ]) 
Y ( ' l) [4 fc,m) ]|GW I Y( £ ) [\i e ' v) ]) (Y [h) [A^' m) ]|GW |Y^ [a£, M) 




gn[X^]d^/kl+al 





(47) 



where the first term corresponds to the h part of the 
operator G and the second one to the e part. Here A 7 
and A„ can be positive and negative. 



If we now consider in addition the limit k y — > 0, we can 
deduce the following properties for the reflection matrix. 
Since = = ©y = for i ^ j we immediately 



have that 



^ ee J(e = o,fc J/ = o) = o, 

£( eft )(£ = 0, k y = 0) = £>' fc )(£ = 0, k v = 0) = 0. (48) 



The only part of the reflection operator which does not 
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vanish is connected with the /i-polarization: 

= 0,k y = 0) = -[e^]- 1 ^ \ (49) 

A rather lengthy calculation also allows us to obtain some 
properties of the previous operator matrix elements. We 
only report the most relevant one for the first Brillouin 
zone, i.e. in the limit aop/n -C 1 

g& h) (Z = 0,k v =0)-loc-a o . (50) 

The solid lines in figure[5]are the numerically computed 
matrix elements of the reflection operator of the metal- 
lic grating corresponding to the zeroth order reflection 
(7 = ^ = 0) for the first seven Matsubara frequencies, 
shown only for k y — 0. The Drude model was used with 
ujp = 8.39 eV, 7 = 0.043 eV, while the grating geom- 
etry is: pi = 160 nm, P2 — 90 nm, and d = 216 nm. 
The behavior of these reflection amplitudes is in agree- 
ment with the predictions made above. The dashed lines 
show the corresponding matrix elements or a flat metal- 
lic surface (Fresnel coefficients), using the same optical 
parameters. From the figure it is clear that the grating 
is less specularly reflecting than a flat surface. At large 
wavevectors the grating reflection amplitudes behave dif- 
ferently with respect to the plane surface: while for the 
plane they asymptotically reach a horizontal line, for the 
grating they have a finite negative slope. However, at 
small wavevectors the behavior of the reflection ampli- 
tudes for the grating and for the flat surface is similar 
(except for the zeroth Matsubara hh reflection ampli- 
tude). This suggest that in this limit the grating may be 
described using an effective medium approximation, in 
which the reflection matrices of the grating are approxi- 
mated by Fresnel coefficients for an homogeneous planar 
interface with an effective permittivity e e ff(w). A fit of 
the numerical results for the grating in figure [5] for the 
h— and e-polarization to Fresnel coefficients gives an ef- 
fective Drude permittivity with a reduction of the plasma 
frequency of about 7.8 times for h polarization and 2.2 
times for the e polarization. The effective dissipation rate 
decreases more for the e-polarization than for the h-one 
(1.4 times against 1.2). 

Finally, let us emphasize that our method allows us 
to deal with the zero Matsubara frequency (6=0 = 0) 
analytically, without resorting to any limiting procedure, 
such as approximating £; = o by a large wavelength mode 
(as used in, for example, pQ). It also avoids problems 
related to the Gibbs phenomenon which complicates the 
calculation, especially for metallic structures. This phe- 
nomenon refers to the oscillations that occur when a 
piecewise discontinuous function, such as our permittiv- 
ity e{x;ui), is approximated by a finite Fourier series. Its 
impact increases with the magnitude of the discontinu- 
ity and, therefore, becomes a more serious issue at low 
frequencies. Instead, the method described in this paper 
deals with such a discontinuity exactly, eliminating de 
facto all problems relate with a Fourier decomposition of 
the permittivity profile. 
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FIG. 5: Matrix elements of the reflection operator correspond- 
ing to the zeroth order reflection. The full lines are the result 
for the grating for the first seven Matsubara frequencies. The 
dashed lines show the equivalent matrix elements for the flat 
surface (Fresnel coefficients), (a): zeroth order reflection co- 
efficients for the /i-polarization. The first seven Matsubara 
frequencies are shown (zeroth to the sixth from the top to 
the bottom), (b): zeroth order reflection coefficients for the 
e-polarization. Once again, the first seven Matsubara frequen- 
cies are shown. The zeroth frequency term is zero, while the 
other ones (first to the sixth from the bottom to the top) de- 
crease in absolute value. The depth of the grating is d — 216 
nm, and the remaining parameters are the same as in fig. [2] 



VII. THE CASIMIR INTERACTION 



In this section we use our quasi-analytical modal ap- 
proach to write down the Casimir free energy ([I]) be- 
tween two vacuum-separated, lamellar gratings facing 
each other. We discuss how to generalize the formal- 
ism to non-lamellar gratings and multilayered periodic 
structures. Finally, we numerically compute the Casimir 
pressure between a flat gold plate parallel to a gold grat- 
ing, and discuss the asymptotic behaviors at large and 
small distances. 
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A. Two lamellar gratings 

Let us consider two vacuum-separated lamellar ID 
gratings. The Casimir pressure between them can be 
obtained by taking the a derivative of the Casimir free 
energy Q, 



show that this matrix is the pivotal matrix 022 mp of the 
composite system formed by the left grating, the vacuum 
region, and the right grating. This is particularly evident 
if one writes it as follows 

©22 mp = 3^ (m)t ' G L {h L ) ■ G (v) (a) ■ G R (h R ) ■ y^ {m \ (55) 



P(a) 



°° />oo r^/p 

V / dk y / da 

1=0 Jo Ja 



x d a log det [l - g L ■ T^ {v) {a) ■ T^ R ■ T^ v) ( 



(51) 



where we have already performed the trace over the spa- 
tial degrees of freedom and used the parity properties of 
the reflection operators 



where 

~(56) 

is the vacuum (propagator) operator. In the limit a — > 
oo (infinitely separated gratings) the second term in the 
numerator of (53) vanishes, and hence we can write (51) 

as 



t L = -^r 1 • e 2 L 1; T^ R = 0f a • [Qg]-\ (52) P(a) dk v da d a log 

P l=0 Jo Jo 



and of V^ v \a). As we discussed above, the calculation 
of the reflection matrices requires the inversion of the 
pivotal matrix, which can be an expensive and not accu- 
rate numerical operation. It is however possible to avoid 
this inversion and derive the Casimir free energy. Indeed, 
using (52) in (51) we can write 



log det [l - £ L • 7>W • -R, R • 7>M 



det 


©2i -t (V \ 


-a) 


©f 2 + ©22 


.%(v)(-a).e R 


dct 


© 2 L 2 




©i 





log 



(53) 

where the different theta-matrices for the left (L) and 
right (R) gratings can be obtained from ([29]), namely 



u 21 



e 



2- 



3>( m >t -G L (d L ) -g} v \ 



0fa = • G R (d R ) ■ % 



m) 



e 



J2 = ^ ■G R (d R )-^ m \ 



(54) 



where dh and d R are the depths of the left and 
right gratings, respectively. The interpretation 
of (53) is particularly simple in terms of modes. 
Indeed, as we discussed above, the determinant 
of the pivotal matrix gives the resonance of the 



system, and hence det 



e 



det [9 



aa] det 



V^X-a) det [9 



7> (»)(-£ 

gives 



nances of the two isolated gratings. 



u 22 



the reso- 
The factor 



det 



represents the contribution 

of the continuum of electromagnetic vacuum modes 
hitting the gratings [30]. Similarly, the determinant of 
©2i ' £ (u) (-a) • ©fa + ©^2 ' T^ v) {-a) ■ ©£, gives the 
coupled modes of the two gratings. Indeed, one can 



^ 00 />oo r^/p 

= — / dky / da Q d a log det 
P l=0 Jo Jo 



det [6 c 2 ° mp (a)] 
dct [6£ mp (a->oo)] 

■y-comp 



(57) 



where we used the definition of the transmission operator 



in terms of the pivotal matrix given in eq.(31 ). 

Before concluding this subsection let us discuss the im- 
pact of mode degeneracy on the Casimir interaction. In 
the case of two degenerate eigenvalues new expressions 
for the eigenfunctions U^ s > must be found using standard 
techniques. Although possible and not mathematically 
involved, this is however of no use for the evaluation of 
the Casimir pressure. Indeed, one can show that this will 
require the modification of the integrand of the previous 
expression in a Lebesgue null measure ensemble of points, 
without changing the final result. 



B. Generalization to non-lamellar gratings 

The simple physical reasoning behind the previous re- 
sults allows us to generalize the calculation to multilay- 
ered periodic structures and non-lamellar gratings, which 
can be approximated, by slicing them, as multilayered pe- 
riodic structures of individual lamellar gratings [311 132] . 
Consider, for example, two non-lamellar gratings facing 
each other and separated by vacuum. The composite 
system is bounded by two homogeneous bulk media ttil 
and m R . The pivotal matrix for the composite system is 
clearly given by 

0™ mp = y(™*)t.(JjGf)-GW.(JjGf).y (m,,) , (58) 



where G i ' is the propagator for the i-th lamellar slice 
of the non-lamellar left or right grating. 
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C. Large distance asymptotic expression 

Despite the complexity of the previous expressions, it 
is possible to derive a close expression for the Casimir 
free energy between gratings in the asymptotic limit of 
large distances, a — !• oo. Since the operator 7^(a) is a 

diagonal matrix with decreasing exponentials e _aA " as 
matrix elements, it follows that, for any fixed Matsubara 
frequency, the eigenvalue with v = is the one that gives 
the slowest decrease as a grows. This value of v corre- 
sponds to the zeroth order of reflection in the standard 
Rayleigh formalism for scattering from periodic struc- 
tures. Therefore, at large distances we can keep only con- 

(v) 

tributions arising from the A„J. eigenvalue, and approxi- 



mate 1 — 
where k 



i- 



' Hi ~ — Iko, 

oo e 



K y 



^]oo-[^ 
q. The subscript 00 indicates 

that only the 7 = ^ = block of the 1Z matrices is con- 
sidered. For the same reason, at distances large enough 
the dominant contributions to the Casimir free energy 
comes from ao rs and k y ss 0. We know already that 
when k y — > the e and h polarizations decouple, which 
implies that the submatrices [^ L ]oo and [7^ fl ]oo become 
diagonal. Hence, in this large distance limit, we approx- 
imate the pressure as 



P(a) 



E 

2=0 



oo 



r/p 
dao 



d a {log [l-[^ L ] e e) [^ 



fil(ee) -2na 



J 00 



log 



1 Lp: Joo Li^ Joo e 



}, (59) 



which is formally equivalent to the integrand of the Lif- 
shitz formula for parallel planes. At large distances, the 
zeroth Matsubara frequency (I — 0) dominates the above 
summation, which implies that P(a) is proportional to 
— fcsTa~ 3 , as in the plane-plane case. The proportion- 
ality factor depends on the value of the reflection am- 
plitudes in the limit k y « and ao ~ 0. For metallic 
gratings described by the Drude model, we have seen 

l(ee) ^-0,fc a =0) =0, 
Therefore 



while [K]ffi(€ = 0,k v = 0,a -> 0) = 1. 
for Drude parallel plates, the prefactor is C(3)/87r. 



above (see {48} and ©) that [TZ] { 00 >(£ 

-+0) 



as 



D. Numerical results 

In this subsection we will focus on the Casimir inter- 
action between a gold lamellar grating parallel to a gold 
flat surface. In principle, our modal approach can treat 
this problem almost fully analytically, requiring numerics 
only for finding the roots of the transcendental equation 
( 34 ) to determine the eigenvalues for the grating region. 



However, from the practical point of view, we are also 
forced to truncate the matrices and the series, to numer- 
ically evaluate integrals, and to deal with convergency 
issues. We address these issues in what follows. 
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FIG. 6: Casimir pressure between a metallic grating and a 
metallic plane, computed using our quasi-analytical modal 
approach. At large separation the pressure tends towards the 
value C{3)k B T/(8na 3 ) (dashed curve). At short separations 
the pressure is oc a~ 3 because of the finite grating conductiv- 
ity. The prefactor used for the dotted curve is the one for the 
plane-plane case multiplied by the filling factor f ~ P2/P (see 
definition of Palling in the text). The Drude parameters are 
ui p — 8.39 eV and 7 = 0.043 eV. The geometrical parameters 
of the grating are: width of the grooves p\ — 160 nm, width of 
the teeth P2 = 90 nm, and height d = 216 nm. Temperature 
is set to T = 300 K. 



The size of the theta- and scattering matrices is set 
by the number of eigenvectors N mayi one keeps to de- 
scribe the fields in the homogeneous regions (equivalent 
to the Rayleigh orders) . This number will be always odd 
because we will truncate the Rayleigh expansion sym- 
metrically with respect to the zeroth order. The corre- 
sponding matrices will be block matrices with dimension 
(2iV max ) x (2A r max ) (the factor 2 comes from the two po- 
larizations) . The expression of the grating operator ( 30 ) 
is formally independent of the truncation order N max , 



and the series defining it could be truncated at a different 
value, say M max . Numerical studies show, however, that 
for the reflection matrix the best convergency is obtained 
when M max = N max . This can be physically understood 
from a argument of dimensionality matching between the 
Hilbert spaces describing the field inside the grating and 
in the homogeneous regions. This is particular clear at 
high frequency where a one-to-one correspondence be- 
tween grating eigenvectors and vacuum eigenvectors is 
required to satisfy the high-frequency transparency. Our 
numerical studies show that, for our choice of optical and 
geometrical parameters (plasma frequency 8.39 eV, dissi- 
pation rate 0.043 eV, p\ — 160 nm, P2 = 90 nm, d = 216 
nm) the first eleven modes (M max = N max = 11) for 
the e- and for /i-polarization are enough for the theta- 
matrices (and, hence, the reflection matrices) to converge 
for all values of £, k y , and ao relevant in the numerics. 
For our configuration, higher modes would correspond 
to values much larger than the plasma frequency, for 
which the metal is almost transparent (see fig|2|. Since 
the magnitude of the eigenvalues decreases with the in- 
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FIG. 7: Plane-plane (dashed) and plane-grating (solid) 
Casimir pressure normalized by -Pmiing with the respective fill- 
ing factor (/ = 1 for the plane-plane and / = 90/250 for the 
plane-grating). The two vertical lines are located at distances 
corresponding to the plasma wave-length (A p = 2-n /ui p ) and 
half of the thermal wave-lenght (At = ftc/2fcsT). The transi- 
tion from the a -3 non-retarded behavior to the a~ 4 retarded 
behavior for the plane-plane configuration occurs much faster 
than for the plane-grating one. The large distance a -3 ther- 
mal regime is however not affected by the grating and starts 
roughly at the same point for both configurations. Parame- 
ters are the same as in previous figures. 



Casimir pressure goes as a~ 3 for a <C A p (non-retarded 
van der Waals regime), as a~ 4 for A p < a < Ay (retarded 
regime), and again as a -3 for a > At (thermal regime), 
where X p — 2ttc/lu p is the plasma wavelength (« 147 nm 
in our case) and At = frc/ksT is the thermal wavelength 
(~ 7 /im at T = 300 K). The behavior at short distances 
can also be interpreted as resulting from the non-retarded 
interaction between surface plasmons [HI EH EO [33J IS] • 
On the other hand, a grating is known for modifying 
the electromagnetic near field, by affecting the behavior 
of surface plasmon modes and effectively increasing the 
plasma wavelength (as we discussed above). Therefore, 
one expects a wider non-retarded regime for the case of 
metallic gratings. Figure [7] shows the grating-plane and 
the plane-plane pressure normalized by Pfiiii ng with the 
respective filling factors (/ = 1 for the plane-plane con- 
figuration). As expected, the transition from the a -3 
to the a -4 behavior happens at a larger distance for the 
plane-grating than for the plane-plane configuration (i.e., 
the grating-plane case has a wider non-retarded regime) . 
The same figure also shows that, on the contrary, the 
thermal regime is not affected and starts roughly at the 
same point for both structures. 



VIII. CONCLUSIONS 



verse of the grating parameters (see (37 1, (38 1, and (|4l[)), 
more (less) modes will be required for gratings with larger 
(smaller) geometrical features. 



The calculation of the Casimir pressure in eq. ( 57 ) also 



requires the evaluation of two integrals over wave vectors. 
The integration is performed using a 30 points Gauss- 
Legendre quadrature scheme for cvo, and a 20 points 
Gauss-Laguerre quadrature scheme for k y . Numerical 
checks show that for the zeroth Matsubara frequency the 
agreement with a Montecarlo calculation is better that 
1 % for 100 nm < a < 5 fim, and better that 3 % for 
5 /Ltm < a < 10 /jm. The agreement greatly improves for 
higher Matsubara frequencies. The Matsubara series was 
truncated at 41 terms. At the distance of a = 50 nm, the 
total result changes by less than 1 % in going from 37 to 
41 Matsubara terms. 

Figure [6] shows the result of the numerical evaluation 
of the Casimir pressure obtained from (57 1. As a check 



of our prediction we also show the large distance asymp- 
totic expression discussed in the previous section (dashed 
line). The dotted line represents the short distance plane- 
plane asymptotic behavior multiplied by the filling factor 
(/ = Pa/p) [9 , Pmiing(a) = -1.79/w p /ic7r/720a 3 . The 
good agreement between the full line and the dotted line 
in Fig. [6] indicates that at short distances the plane- 
grating Casimir pressure is substantially less than the 
plane-plane pressure mainly due to geometrical effects. 

Finally, we briefly address the influence of finite con- 
ductivity and temperature for Casimir interactions in- 
volving gratings. In the plane-plane configuration the 



In summary, we have developed a quasi-analytical 
modal approach to computing Casimir interactions in- 
volving ID lamellar gratings. The method can be gener- 
alized to more complex nanostructures by approximating 
them via slicing as a collection of multilayered lamellar 
gratings [3TJ [35] . The key features of our method is that 
the eigenmodes of the grating can be solved for analyti- 
cally, while the cigcnfrcqucncics arc solutions to a simple 
transcendental equation (34). Apart from these funda- 



mental aspects, we have also presented an approach to 
calculate the Casimir interaction without resorting to any 
matrix inversion that avoids several potential numerical 
instabilities, improves the precision of the numerical re- 
sults, and can be used in other non-modal frameworks. 
We studied analytically the form of the eigenvalues in 
some specific limiting cases, and discussed their impact 
on the scattering operators and on the Casimir interac- 
tion. By analyzing the mode structure at real frequen- 
cies, this formalism can also be applied to study other 
fluctuation-induced interactions (thermal emission, near- 
field heat transfer, etc.). 
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